setwd("/users/brucedesmarais/Dropbox/SenPE/Data")
 library(foreign)
 library(lme4)
 library(MASS)
 sennbreg <- read.dta("sennbreg.dta")
 fowPlus <- read.dta("fowler_plus.dta")
 
 fowId <- paste(fowPlus$labels,fowPlus$congress,sep="")
 senId <- paste(sennbreg$labels,sennbreg$congress,sep="")
 senfow <- match(senId,fowId)
 
sennbreg[,c("nc","ncl","dc","dcl")] <- fowPlus[senfow,c("nc","ncl","dc","dcl")]
 
sennbreg$id <- 1:nrow(sennbreg)
chairs <- read.csv("Chairs.csv",stringsAsFactors=F)

last_name <- NULL
for(i in 1:nrow(sennbreg)){
	namei <- sennbreg$labels[i]
	lasti <- toupper(strsplit(namei,split=",")[[1]][1])
	last_name <- c(last_name,lasti)
	}

chair <- numeric(length(last_name))

congs <- 97:105
for(t in congs){
	chairst <- subset(chairs,chairs$congress==t)
	chairst <- c(chairst$chair1,chairst$chair2)
	chairst <- chairst[which(nchar(chairst)>1)]
	for(i in 1:length(chairst)){
		chairi <- chairst[i]
		chairi <- toupper(strsplit(chairi,split=",")[[1]][1])
		ind <- which(last_name==chairi & sennbreg$congress==t)
		chair[ind] <- 1
		}
	}

# Most Central Senator Table

CentTab <- NULL
for(c in congs){
	centc <- subset(sennbreg,sennbreg$congress==c)$nc
	idLabel <- subset(sennbreg,sennbreg$congress==c)$labels
	ord <- order(-centc)
	top <- idLabel[ord[1:5]]
	CentTab <- cbind(CentTab,top)
	}

colnames(CentTab) <- congs

CentTab <- CentTab[,c(5,7,9)]

library(xtable)
xtable(CentTab)

sennbreg$chair <- chair

## Pooling ##

ideol_extrm <- numeric(nrow(sennbreg))
majority <- numeric(nrow(sennbreg))
ucong <- unique(sennbreg$congress)
for(i in 1:length(ucong)){
	congi <- which(sennbreg$congress==ucong[i])
	ideoli <- sennbreg$ideol1[congi]
	partyi <- sennbreg$party[congi]
	extri <- abs(ideoli-median(ideoli))
	majority[congi] <- as.numeric(partyi==median(partyi))
	ideol_extrm[congi] <- extri
	}

fowler <- read.csv("SH.csv")
	
sennbreg$ideol_extrm <- ideol_extrm
sennbreg$obs <- 1:nrow(sennbreg)
sennbreg$majority <- majority
sennbreg$l_connectedness[which(is.na(sennbreg$l_connectedness))] <- 0

fow96ind <- which(sennbreg$congress==97)
fowInlast <- match(paste(sennbreg$congress,"_",sennbreg$ids,sep=""),paste(fowler$congress-1,"_",fowler$ids,sep=""))
sennbreg$pb_l <- fowler$pb[fowInlast]
sennbreg$pb_l[which(is.na(sennbreg$pb_l))] <- 0
sennbreg$pa_l <- fowler$pa[fowInlast]
sennbreg$pa_l[which(is.na(sennbreg$pa_l))] <- 0
sennbreg$sponsored_l <- fowler$sponsored[fowInlast]
sennbreg$sponsored_l[which(is.na(sennbreg$sponsored_l))] <- 0

library(lme4)

LRtest <- function(mod2,mod1){
	## Mod2 has more df
	ll2 <- logLik(mod2)
	df2 <- attr(ll2,"df")
	ll1 <- logLik(mod1)
	df1 <- attr(ll1,"df")
	tstat <- 2*(as.numeric(ll2)-as.numeric(ll1))
	dftest <- df2-df1
	pval <- 1-pchisq(tstat,dftest)
	list(test=tstat,df=dftest,pval=pval) 
} 

### Models of Bill Passage

# Baseline
baseline <- lmer(PassS~ideol_extrm+chair+log(Sponsored+.0001)+seniority*majority+(1|obs)+as.factor(congress)+(1|ids), data=sennbreg,subset=!is.na(sennbreg$connectedness) & !is.na(sennbreg$l_connectedness),family=poisson,REML=F)

# Cosponsorship
cosponsor <- lmer(PassS~connectedness+ideol_extrm+chair+log(Sponsored+.0001)+seniority*majority+(1|obs)+as.factor(congress)+(1|ids), data=sennbreg,subset=!is.na(sennbreg$connectedness) & !is.na(sennbreg$l_connectedness),family=poisson,REML=F)

# Both Additive
additive <- lmer(PassS~nc+connectedness+ideol_extrm+chair+log(Sponsored+.0001)+seniority*majority+(1|obs)+as.factor(congress)+(1|ids), data=sennbreg,subset=!is.na(sennbreg$connectedness) & !is.na(sennbreg$l_connectedness),family=poisson,REML=F)

# Interaction
both <- lmer(PassS~nc*ncl+connectedness+ideol_extrm+chair+log(Sponsored+.0001)+seniority*majority+(1|obs)+as.factor(congress)+(1|ids), data=sennbreg,subset=!is.na(sennbreg$connectedness) & !is.na(sennbreg$l_connectedness),family=poisson,REML=F)

# both interactions
interact <- lmer(PassS~nc*ncl+connectedness*l_connectedness+ideol_extrm+chair+log(Sponsored+.0001)+seniority*majority+(1|obs)+as.factor(congress)+(1|ids), data=sennbreg,subset=!is.na(sennbreg$connectedness) & !is.na(sennbreg$l_connectedness),family=poisson,REML=F)

models <- list(baseline,cosponsor,additive,both, interact)
#### Reporting from Committee

# Baseline
baselineC <- lmer(ReportS~ideol_extrm+chair+log(Sponsored+.0001)+seniority*majority+(1|obs)+as.factor(congress)+(1|ids), data=sennbreg,subset=!is.na(sennbreg$connectedness) & !is.na(sennbreg$l_connectedness),family=poisson,REML=F)

# Cosponsorship
cosponsorC <- lmer(ReportS~connectedness+ideol_extrm+chair+log(Sponsored+.0001)+seniority*majority+(1|obs)+as.factor(congress)+(1|ids), data=sennbreg,subset=!is.na(sennbreg$connectedness) & !is.na(sennbreg$l_connectedness),family=poisson,REML=F)

# Both Additive
additiveC <- lmer(ReportS~nc+connectedness+ideol_extrm+chair+log(Sponsored+.0001)+seniority*majority+(1|obs)+as.factor(congress)+(1|ids), data=sennbreg,subset=!is.na(sennbreg$connectedness) & !is.na(sennbreg$l_connectedness),family=poisson,REML=F)

# Interaction
bothC <- lmer(ReportS~nc*ncl+connectedness+ideol_extrm+chair+log(Sponsored+.0001)+seniority*majority+(1|obs)+as.factor(congress)+(1|ids), data=sennbreg,subset=!is.na(sennbreg$connectedness) & !is.na(sennbreg$l_connectedness),family=poisson,REML=F)

# both interactions
interactC <- lmer(ReportS~nc*ncl+connectedness*l_connectedness+ideol_extrm+chair+log(Sponsored+.0001)+seniority*majority+(1|obs)+as.factor(congress)+(1|ids), data=sennbreg,subset=!is.na(sennbreg$connectedness) & !is.na(sennbreg$l_connectedness),family=poisson,REML=F)

modelsC <- list(baselineC,cosponsorC,additiveC,bothC, interactC)



set.seed(5)
sim_both <- simulate(both,5000)
quants <- seq(0.1,0.99,length=500)
y <- attributes(both)$y
qobs <- quantile(y,prob=quants)
qsim <- quantile(c(as.matrix(sim_both)),quants)

# qq-plot
# histograms with all means overlayed

# qq-plot
pdf("/users/brucedesmarais/Dropbox/SenPE/Tex/qq_fit.pdf",height=3.5, width=4, family="Times", pointsize=13.5)
par(las=1,mar=c(4,5,.5,.5),cex.lab=1.25)
plot(qsim,qobs,ylim=c(0,15),xlim=c(0,15),col="white", ylab="Observed Quantile",xlab="")
abline(0,1,lwd=2,col="grey50")
title(line=2.25,xlab="Simulated Quantile")
points(qobs,qsim,pch=4)
dev.off()

group_means <- function(object,sims, grp){
	group <- attributes(object)$flist[,grp]
	y <- attributes(object)$y
	ugroup <- unique(group)
	grp_mean <- numeric(length(ugroup))
	sim_mean <- matrix(0,length(ugroup),ncol(sims))
	for(i in 1:length(ugroup)){
		whichi <- which(group==ugroup[i])
		grp_mean[i] <- sum(y[whichi])
		for(j in 1:ncol(sims)){
			sim_mean[i,j] <- sum(sims[whichi,j])
			}
		}
		list(obs=cbind(ugroup,grp_mean),sim=sim_mean)
	}
	
group_means2 <- function(object,sims, grp){
	group <- grp
	y <- attributes(object)$y
	ugroup <- unique(group)
	grp_mean <- numeric(length(ugroup))
	sim_mean <- matrix(0,length(ugroup),ncol(sims))
	for(i in 1:length(ugroup)){
		whichi <- which(group==ugroup[i])
		grp_mean[i] <- sum(y[whichi])
		for(j in 1:ncol(sims)){
			sim_mean[i,j] <- sum(sims[whichi,j])
			}
		}
		list(obs=cbind(ugroup,grp_mean),sim=sim_mean)
	}	

## Histogram illustrating overall mean
pdf("/users/brucedesmarais/Dropbox/SenPE/Tex/hist_overall.pdf",height=3.5, width=4, family="Times", pointsize=13.5)
par(las=1,mar=c(4,5,.5,.5),cex.lab=1.25)
hist(apply(sim_both,2,mean),freq=F,col="grey75",border="grey55",xlab="",main="")
abline(v=mean(y),lwd=2.5)
title(line=2.25,xlab="Mean Bills Passed")
dev.off()


sens_means <- group_means(both,sim_both,2)
senMu <- sens_means[[1]][,2]
senSim <- t(sens_means[[2]])
axis(1,line=1)
rug(senMu,ticksize=.05,line=1,col=rgb(0,0,0,75,maxColorValue=255))

## Histogram illustrating Senator Total
pdf("/users/brucedesmarais/Dropbox/SenPE/Tex/hist_senator.pdf",height=3.5, width=4, family="Times", pointsize=13.5)
par(las=1,mar=c(5,5,.5,.5),cex.lab=1.25)
hist(senSim,xlim=c(0,400),breaks=50,freq=F,col="grey75",border="grey55",xlab="",main="",xaxt="n",ylab="")
axis(1,line=1)
rug(senMu,ticksize=.1,line=1,col=rgb(0,0,0,75,maxColorValue=255))
title(line=3.25,xlab="Total Bills/Senator")
title(line=3.75,ylab="Density")
dev.off()

cong_means <- group_means2(both,sim_both,sennbreg[!is.na(sennbreg$connectedness) & !is.na(sennbreg$l_connectedness),"congress"])
congMu <- cong_means[[1]][,2]
congSim <- t(cong_means[[2]])

## Histogram illustrating Congress Total
pdf("/users/brucedesmarais/Dropbox/SenPE/Tex/hist_cong.pdf",height=3.5, width=4, family="Times", pointsize=13.5)
par(las=1,mar=c(5,6,.5,.5),cex.lab=1.25)
hist(congSim,breaks=40,freq=F,col="grey75",ylab="",border="grey55",xlab="",main="",xaxt="n")
axis(1,line=1)
rug(congMu,ticksize=.1,line=1,col=rgb(0,0,0,150,maxColorValue=255))
title(line=3.25,xlab="Total Bills/Congress")
title(line=4,ylab="Density")
dev.off()

### Interpretation of the Main Effect of Interest ###
library(mvtnorm)
set.seed(5)
sds <- sqrt(unlist(VarCorr(both)))
re <- 0
for(i in 1:length(sds)){
	re <- rnorm(1000,sd=sds[i])
	}

nc <- 2
ncl <- 3
nc_ncl <- 18

X <- attributes(both)$X
nclDat <- X[,ncl]
qncl <- quantile(nclDat,prob=c(.125,.375,.625,.875))

quarts <- cbind(quantile(nclDat,prob=c(0,.25,.5,.75)),quantile(nclDat,prob=c(.25,.5,.75,1)))

efflist <- list()
for(i in 1:length(qncl)){
X <- attributes(both)$X
X <- subset(X,X[,ncl] >= quarts[i,1] & X[,ncl] <= quarts[i,2])
Xoth <- X[,-c(nc,ncl,nc_ncl)]
beta <- attributes(both)$fixef
betaOth <- beta[-c(nc,ncl,nc_ncl)]
Xmu <- apply(Xoth,2,median)

lpOth <- t(betaOth)%*%cbind(Xmu)

beta_mu <- beta[c(nc,ncl,nc_ncl)]
beta_cov <- vcov(both)[c(nc,ncl,nc_ncl),c(nc,ncl,nc_ncl)]
beta_cov <- round(as.matrix(beta_cov),9)

betas <- rmvnorm(10000,beta_mu,sigma=beta_cov)

nclDat <- X[,ncl]
qncl <- quantile(nclDat,prob=c(.25,.5,.75,1))
ncDat <- X[,nc]
xseq <- seq(min(ncDat),max(ncDat),length=100)

set.seed(13453453)
	ncli <- qncl[i]
	effmati <- matrix(0,10000,100)
	for(j in 1:length(xseq)){
		lpij <- lpOth
		ncj <- xseq[j]
		for(k in 1:nrow(betas)){
			lpijk <- lpij+ncj*betas[k,1]+ncli*betas[k,2]+ncj*ncli*betas[k,3]+re
			lambda <- exp(lpijk)
			#yijk <- rpois(10000,lambda)
			#muijk <- mean(yijk)
			effmati[k,j] <- mean(lambda)
			}
		print(t(c(i,j)))
		}
	efflist[[i]] <- effmati

filei <- paste("/users/brucedesmarais/Dropbox/SenPE/Tex/effect",i,".pdf",sep="")
pdf(filei,height=3.5, width=4, family="Times", pointsize=13.5)
par(las=1,mar=c(4,4,.5,.5),cex.lab=1.25)
plot(xseq,apply(efflist[[i]],2,mean),type="l",ylim=c(0,7),ylab="Mean Bills Passed",xlab="")
polygon(c(xseq,rev(xseq)),c(apply(efflist[[i]],2,quantile,prob=0.975),rev(apply(efflist[[i]],2,quantile,prob=0.025))),col="grey65",border="grey65")
lines(xseq,apply(efflist[[i]],2,mean),lwd=2.5)
rug(ncDat,ticksize=0.075,col=rgb(0,0,0,25,maxColorValue=255))
title(line=2.25,xlab="Press Event Centrality")
dev.off()
	}


### Make table, baseline, cosponsor_base, network_base, both
# VarCorr-random effects
models <- list(baseline,cosponsor,additive,both, interact)
effnames <- NULL
for(i in 1:length(models)){
	effnames <- c(effnames,rownames(attributes(summary(models[[i]]))$coefs))
	}
ueffs <- sort(unique(effnames))

effMat <- matrix(0,length(ueffs)*2,5)
reMat <- matrix(0,2,5)
fitMat <- matrix(0,3,5)

sig <- qnorm(0.975)
for(i in 1:length(models)){
	effi <- attributes(summary(models[[i]]))$coefs
	coefi <- effi[,1]
	sei <- effi[,2]
	sigi <- which(abs(coefi/sei)>sig)
	coefi <- round(coefi,3)
	coefi[sigi] <- paste(coefi[sigi],"*",sep="")
	sei <- paste("(",round(sei,3),")",sep="")
	effnamesi <- rownames(effi)
	indi <- match(effnamesi,ueffs)
	effMat[indi*2-1,i] <- coefi
	effMat[indi*2,i] <- sei
	reMat[,i] <- unlist(VarCorr(models[[i]]))
	fitMat[,i] <- c(logLik(models[[i]]),AIC(models[[i]]),LRtest(models[[i]],models[[1]])$pval)
	}
rnEm <-1:(length(ueffs)*2)
rnEm[(1:length(ueffs))*2-1] <- ueffs
rnEm[which(is.na(rnEm))] <- (1:22)[which(is.na(rnEm))]

rownames(effMat) <- rnEm
rownames(reMat) <- c("Observation","Senator")
rownames(fitMat) <- c("LL","AIC","BIC")
nMat <- matrix(rep(c(903,189,9),5),3,5)
rownames(nMat) <- c("obs","sen","cong")

allMat <- rbind(effMat,round(reMat,3),round(fitMat,3),nMat)
library(xtable)
xtable(allMat)	

### Make table, baseline, cosponsor_base, network_base, both
# VarCorr-random effects
models <- list(baseline,cosponsor,additive,both, interact)
effnames <- NULL
for(i in 1:length(models)){
	effnames <- c(effnames,rownames(attributes(summary(models[[i]]))$coefs))
	}
ueffs <- sort(unique(effnames))

effMat <- matrix(0,length(ueffs)*2,5)
reMat <- matrix(0,2,5)
fitMat <- matrix(0,3,5)

sig <- qnorm(0.975)
for(i in 1:length(models)){
	effi <- attributes(summary(models[[i]]))$coefs
	coefi <- effi[,1]
	sei <- effi[,2]
	sigi <- which(abs(coefi/sei)>sig)
	coefi <- round(coefi,3)
	coefi[sigi] <- paste(coefi[sigi],"*",sep="")
	sei <- paste("(",round(sei,3),")",sep="")
	effnamesi <- rownames(effi)
	indi <- match(effnamesi,ueffs)
	effMat[indi*2-1,i] <- coefi
	effMat[indi*2,i] <- sei
	reMat[,i] <- unlist(VarCorr(models[[i]]))
	fitMat[,i] <- c(logLik(models[[i]]),AIC(models[[i]]),LRtest(models[[i]],models[[1]])$pval)
	}
rnEm <-1:(length(ueffs)*2)
rnEm[(1:length(ueffs))*2-1] <- ueffs
rnEm[which(is.na(rnEm))] <- (1:22)[which(is.na(rnEm))]

rownames(effMat) <- rnEm
rownames(reMat) <- c("Observation","Senator")
rownames(fitMat) <- c("LL","AIC","BIC")
nMat <- matrix(rep(c(903,189,9),5),3,5)
rownames(nMat) <- c("obs","sen","cong")

allMat <- rbind(effMat,round(reMat,3),round(fitMat,3),nMat)
library(xtable)
xtable(allMat)	

### Make table, baseline, cosponsor_base, network_base, both
# VarCorr-random effects
models <- list(baselineC,cosponsorC,additiveC,bothC, interactC)
effnames <- NULL
for(i in 1:length(models)){
	effnames <- c(effnames,rownames(attributes(summary(models[[i]]))$coefs))
	}
ueffs <- sort(unique(effnames))

effMat <- matrix(0,length(ueffs)*2,5)
reMat <- matrix(0,2,5)
fitMat <- matrix(0,3,5)

sig <- qnorm(0.975)
for(i in 1:length(models)){
	effi <- attributes(summary(models[[i]]))$coefs
	coefi <- effi[,1]
	sei <- effi[,2]
	sigi <- which(abs(coefi/sei)>sig)
	coefi <- round(coefi,3)
	coefi[sigi] <- paste(coefi[sigi],"*",sep="")
	sei <- paste("(",round(sei,3),")",sep="")
	effnamesi <- rownames(effi)
	indi <- match(effnamesi,ueffs)
	effMat[indi*2-1,i] <- coefi
	effMat[indi*2,i] <- sei
	reMat[,i] <- unlist(VarCorr(models[[i]]))
	fitMat[,i] <- c(logLik(models[[i]]),AIC(models[[i]]),LRtest(models[[i]],models[[1]])$pval)
	}
rnEm <-1:(length(ueffs)*2)
rnEm[(1:length(ueffs))*2-1] <- ueffs
rnEm[which(is.na(rnEm))] <- (1:22)[which(is.na(rnEm))]

rownames(effMat) <- rnEm
rownames(reMat) <- c("Observation","Senator")
rownames(fitMat) <- c("LL","AIC","BIC")
nMat <- matrix(rep(c(903,189,9),5),3,5)
rownames(nMat) <- c("obs","sen","cong")

allMat <- rbind(effMat,round(reMat,3),round(fitMat,3),nMat)
library(xtable)
xtable(allMat)
